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Quantum Simulation of Phylogenetic Trees 
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Quantum simulations constructing probability tensors of biological multi-taxa in phylogenetic trees 
are proposed, in terms of positive trace preserving maps, describing evolving systems of quantum 
walks with multiple walkers. Basic phylogenetic models applying on trees of various topologies 
are simulated following appropriate decoherent quantum circuits. Quantum simulations of statis- 
tical inference for aligned sequences of biological characters are provided in terms of a quantum 
pruning map operating on likelihood operator observables, utilizing state-observable duality and 
measurement theory. 



Introduction: In the last two decades quantum me- 
chanics has found itself in a situation that could be 
characterized as an epistemological exodus. It has 
expanded its scope and applicability into new fields, 
such as information theory, the theory of computa- 
tion, and even biology, and has addressed fundamental 
problems and procedures of these fields, by means of 
its physical-mathematical conceptual and computational 



apparatu 
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What were previously accepted as quan- 



tum paradoxes and oddities, like quantum entanglement, 
have turned out to be the keys to constructing novel 
computational and communicational algorithms, provid- 
ing the means for launching a new quantum technology. 
In this vein, this paper puts forward a novel application 
of the discipline of quantum computation-information to 
the field of evolutionary phylogenetics^i^. Phylogenetics' 
main task is to construct ancestral relationships (phylo- 
genies), inferred by analyzing statistical data, collected 
for various (morphological or genotypic) kinds of char- 
acters or traits, possessed by selected groups of biolog- 
ical organisms (taxa). This amounts to construction of 
phylogenetic trees with appropriate branching patterns 
and evolutionary lengths, that successfully reproduce 
statistical trends of alignments of sequences of certain 
characters'^'-. Various evolutionary models that compete 
by adjusting their tree vertex transition probabilities, to 
accomplish this computationally NP-hard task^, are then 
assessed by some statistical estimation such as maximum 
likelihood measure^. 

In this work a quantum simulation of phylogenetic 
evolution and inference, is introduced in terms of trace 
preserving maps operating on quantum density matri- 
ces. Basic multi-parametric evolutionary models are sim- 
ulated, and an association between phylogenetic trees 
and quantum circuits is established. Specifically, group- 
based models are associated to quantum walks (QW), 
and the Felsenstein model is shown to be related to post- 
measurement state maps. Finally a quantum simulation 
of the iterative pruning process for estimating maximum 
likelihood of phylogenetic trees, is established in terms 
of quantum measurements of likelihood operator valued 



measures (observables) . 

Notation: Let the character set be S = {0,1,...,A^ — 
1} = {0} U E*. Here is considered to be the "null" 
or no character symbol. Introduce the Hilbert space of 
character states H w hi'^) — span{\i);i G E), of di- 
mension dimH = |S| , and consider the space Lin{H) 
of linear operators acting on H. Examples are the com- 
plete set of projectors Pi = \i) {i\ , i & E, the shift oper- 
ator h \i) = \i +N 1) , with +N addition modulo N, (so 
that = 1), and the space of density matrices 'D{II) C 
Lin{H). A classical (discrete) probability distribution is 
represented as a vector {po,pi,p2, ...pN-i), and the corre- 
sponding quantum stochastic system is represented by a 
diagonal density matrix p = X^ies* Pi^i ^ ^(^); fo'^ bio- 
logical applications we will always assume po = (so that 
in practice the sum runs only over characters i G E*). 
On bipartite systems, the unitary control-not operator 
Ucn G Lin{H ® H) defined as Ucn — X^fces Pk®h'', acts 



asi Uc 



-N ■ 

Splitting, cladogenesis, speciation: The splitting 
operation^ A for given 1-taxon matrix p = Xies'^**-^*' 
is implemented by the adjoint action of Ucn 



(1) 



where pij = PiSij, so Ap is identified with a two-taxon 
density matrix. The control-not gate embedded in vari- 
ous positions in s-fold products of character spaces, e.g. 
-^(Sk-i ^ jj^^ ^ -^^s-k-i^ provides the way to construct 

s-taxon phylogenetic trees of various topologies^. 

Phyletic evolution, anagenesis: For an s-taxon density 
matrix p = J2i^ i^es* Pti-.-.t^Ph <8) . . . <8) -Pi, , a suitable 
local unitary U = (^^^i Ui G Lin{H)^^, formalizes the 
phyletic evolution of taxa, when its action is composed 
with the s-fold product of the local diagonalizing map 
f J", where Sd{-) = Xfees Pk{-)Pk, is the completely pos- 
itive trace preserving (CPTP) map that projects out the 
diagonal part of a matrix^, that is a decoherent map. 
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FIG. 1. A 4-t£ixa tree and its simulating quantum circuit. 



Thus we have p^p = £f'{UpU''), where 

P= P^l■■■■i,P^^(S) ■.■(S)P^,, (2) 

ii,...,is6S* 

Pii....is^ Y Pji....jAMi Ms)i,j,-...-i^j^. (3) 

Abbreviating the adjoint action of an operator as 
AdS{-) EE S{-)S\ we say that the map £f'{AdU) thus 
induces a general doubly-stochastic transformation in the 
probability tensor. The Hadamard or entry-wise product 
of matrices defined as {Ao B)ij = AijBij, has been used, 
to obtain- the Markov matrices Mi — UiO U*, which will 
drive evolution on edges of a model phylogenetic tree. 
Below, we make particular choices of U to reflect differ- 
ent types of phylogenetic models. Fig. [T] summarizes the 
preceding discussion by showing a four taxon tree and its 
simulating quantum circuits. 

Phyletic evolution and quantum walks: It has long 
been appreciated that faithful modeling of trait evolu- 
tion in phylogenetics is problematic. As has been re- 
marked, "...Brownian motion is a poor model, and so 
is Ornstein-Uhlenbeck"— . We here present a novel pro- 
posal for the stochastic phyletic evolution of traits via 
quantum simulation employing QWs (see^), operating 
locally on density matrices along edges of trees. This 
is set up as follows. Introduce in additional to char- 
acter Hilbert space H (the "walker" space), at each 
node of phylogenetic tree an auxiliary "coin" Hilbert 
space He « '2(C) = span{\+) ,\—)), and projectors 
P± € Lin{Hc). Evolution now proceeds on joint "walker" 
and "coin" states pc® p via a standard QW conditional 
unitary operator V — {P+ h + P- ^ h^)U €5 1, act- 
ing from He (S) H to itself. One "step" of such a QW 
is realized by the map on the "walker" density matrix, 
viz. p — )■ Syk (p) := TrcV'ipc (g> p)^^^, followed by diag- 
onalization with £d- For s taxa, Eyk = {£f^ o 
For example for the two-taxon case, with k = 2 and 
coin initially in a pure state pc = |c) (c| with |c) ~ |+) 
or |-), we obtain Ey^ip) = J2m7iPmnPm fXi P„, with 

components p,„„ = Y.abP'm-a,n~bqfW^\ where qf^ 
Mj_a-'yMj^a,c > 0, is a probability distribution 



(that is, qi'^^ > 0,J2s1^'^^ ~ 1)' determined by the coin 
tossing unitary U via the Hadamard product M = UoU*. 
The tensor p so obtained, and its multi-taxa generaliza- 
tions, are objects of quantum simulations. Also the diag- 
onalizing map £d can be cast in the form of a CPTP map, 

i.e. £d{p) = Z^fees* PkpPk Z^/tes* QkUkpUl with each 
Qk = thanks to the non- uniqueness of the operator 

sum representation, with unitaries Uk related to projec- 
tors by discrete Fourier transform, Uk — J^i'^^ Pi ''^^'^ 
LJ = exp(j27r/ Below, similar quantum prescriptions 
will be given to the structural maps of standard evolu- 
tionary models. 

Phylogenetic evolutionary models and quantum maps: 
Next we exploit the above considerations in specific cases 
of standard phylogenetic models, namely the so-called 
group-based models (see references^): Jukes-Cantor 
(JC), Kimura two- parameter (K2), Kimura three- 
parameter (K3), and the binary symmetric model (B), 
as well as the Felsenstein model (F)i^. Firstly we give in 
each case a direct Kraus representation of the quantum 
map Er = £do£r- This is followed by a Q W formulation 
using, as above, an additional ancillary "coin" space. 
Let X, Z denote the usual single qubit not and phase 
gates (the Pauli matrices ax, Oz respectively) and Uki — 
(g) X\ for k,l = 0,1. The following propositions 
are verified by direct calculation for operators in ^2(2^*) 
acting on p = Emss* PmPm- 

Proposition K: Let \Y.*\ = 4 and r G {Ki,K2, JC}. We 
have 

i?.(p)-^A(;)C/fei(p){/t^= ^ {MrP)mPm, 



k.l 



Mr{a,b,c)=J2KiUkioU, 



kl 



kl - X! '^kl^^" 
kl 



<XK (4) 



The weights AJ; and corresponding model Markov ma- 
trices Mt are defined as follows. For generic parameters 
define the weights Xki{a,, b, c) as Aqo = 1—a—b—c, Aio = a, 
Aoi = b, All = c, and take the corresponding convex sum 
Af(a,6,c). Then a[^;^^ = Xki{a,b,c), M^k = M{a,b,c), 



A 



{2K) 
kl 



= Xki{a,b,b), M2K — M{a,b,b), and finally 



^ki^^ ^ Xki{a,a,a), Mjc = M{a,a,a). 



□ 



Proposition K' : The CPTP map E^- has, in addition to 
the operator sum representation above, also a QW like 
representation Et-{p) = TvcVripc ® p)^r^, m terms of a 
unitary dilation Vr = (Zfc; Pk® Pi® Uki)Ur ® 1 which 
acts on a composite coin- walker space H(.®H, with four- 
dimensional ancillary space. Here Vr is a control-control- 
Uki operator. For a coin density matrix with spectral 
decomposition pc — Zfe Mfc |cfc) (cfe| , the coin-tossing 

with 



unitary Ur should satisfy {kl\ Ur o U* |c) 
|c) = Za: Pk \ck) a stochastic vector. Also Uki — e 
where Tiki = ^'''[-(^ + 1)1 (g) 1 + kX ® 1 + 11 ^ X]. □ 



kl ' 



Proposition B: Let \Yi* 
Eb{p) = (1 - a)p -f 



— 2. The map £d ° Eb, where 
aX pX\ simulates the binary 
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symmetric model Mb {a) = (1 — a)l+aX acting as 

P = Ernes* PmPm EmeS* {MBp)mPm- □ 

Proposition B' : The "control flip" map Eb is unitarized 
in composite coin-walker space with a two-dimensional 
ancillary space as, Eb{p) — TrcVB{pc ® pWb^ 
with the starting coin state pc = |1) and 
= Val®l-h yT^F(8)X, and y = ZX. □ 

Remark: In the QW picture, the weight parameters A^^"* 
determine non-uniquely, via the unistochastioi^ matrix 
Ut o U*, the coin-tossing matrix Uc, which in turn de- 
termines the [/-quantization of the underlying classical 
walk^^ with evolution matrix Vd = Efe/ ^ki- 

For the Felsenstein model (F)^, quantum simulation 
requires the following ingredients. The model's sta- 
tionary distribution (7ri,7r2, 713,714), ^ tTj = 1, is to be 
used to introduce the observable Itt := 4^.7riPi, with 
Kraus operators Fij = ^/Wj\i) {j\ , i,j £ T, = {1,2,3,4} 

obeying the resolution relation J^ij ^ij^v = jItt- Again 
let p = Ernes* PmPm- By direct calculation we obtain: 

Proposition F: The quantum map implementing the 

Felsenstein model p Ef{p) = Emes* (^^^^ P)™^™ 
given by 

Ep{p) = (1 - a)— ^ F,,pFl + ap, (5) 

where p^, = Tr{J2i,j F^jFijP) = Tr{jlT,p) is a normal- 
ization constant, and the model's stochastic matrix is 
obtained as Alp = (1 — a) ^ ■ j Fij o F^j + al. 

In the framework of quantum measurement theory, 
simulation of the Felsenstein model is interpreted as fol- 
lows. There are two observables: 1^ as above, and 
also if defined analogously in terms of the comple- 
mentary probability distribution {irf ,Trf ,Trf), with 
1 = +1^" forming a non-orthogonal decomposition of 
unity. These observables are measured by means of the 
so called instruments^, which are the two families of 
Kraus generators: the {Fij}} as above, and the anal- 
ogous {F^}lj^Q defined in terms of tt'^ rather than tt 
(see e.g^). The measurement probabilities of the ob- 
servables 1^ and if in the system are pf — Tr{lfp), 
and Pti = Tr^l^rp), and the action of quantum map Ep 
on the density matrix Ep{p) gives the post-measurement 
density matrix for a non-efficient quantum measurement 
for observable Ijr of finite strength^^ . The complemen- 
tary measurement of if is not used. In the uniform limit 
TTj = i then = 1, 1^ = and pj,- — 1, and the model 
reduces to the JC model. 

Quantum estimation of likelihood: Our general frame- 
work also encompasses the quantum estimation of model- 
based tree likelihoods (F)i^, whose numerical calculation 
and optimization provides a major tool for phylogenetic 
inference (for computational heuristics see e.gJ^). Like- 
lihood evaluation has been demonstrated to be a com- 
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FIG. 2. Circuit for pruning map of likelihood operators. 



putationally NP-hard problem^, and it is therefore de- 
sirable to put forward a quantum simulation equivalent. 
In the usual formulation (F) , likelihood vectors are ini- 
tialized at the pendant nodes (leaves) of a tree, and are 
then computed recursively back to the root node, the 
final result being a scalar quantity, the tree likelihood. 
The key operation is that of pruning, that is, of arriv- 
ing at the likelihood for a parent node, say A, by com- 
bining a pair of daughter likelihoods, say B, C, from 
nodes which root two sub-trees. Explicitly, likelihoods for 
daughter nodes B, C are combined to give the parent like- 
lihood Lf = (E,es. (ts)if )(Efees- Mg{tc)L<^), 
where M^-'-^ {tB,c) are stochastic matrices depending on 
branch lengths ts^c specified by the evolutionary model 
employed. Next, an alignment of s taxa over A sites 
is considered. If the characters at site I of the align- 
ment are 4'''4' then likelihoods for the tips of the 
tree (leaf nodes) are initialized to L^'' = 5{k^if'). The 
pruning map is applied recursively at all cherries, and 
then higher up the tree, to arrive at the total tree like- 
lihood L^^} = (i[r''fe)fe=i' which is finally averaged over 
the assumed stationary distribution (7ri) of the model to 

obtain site Vs likelihood L*^'^ = Efe'^'i^tr fe- Foi' the en- 
tire alignment, the tree (log) likelihood is then L(T;w*) 
= maxu, logll^j^ where T denotes the tree topology 
and w* the optimal model (weight) parameters. 

In the quantum simulation introduced here, likeli- 
hoods are regarded as quantum observables, that is op- 
erators in Lin{H), dual to density operators under the 
trace inner product (see above). The likelihood opera- 
tor at node A has components Lf = L^{t\i) — P{i\t), 
where ¥{i\t) is the conditional probability of character 
i G S*, for parameters t = {T,w). Here A — l,2,...,s 
are leaf nodes and A = s + 1, ...,2s — 2, internal (an- 
cestral) nodes. Consider parent and daughter nodes A, 
B and C, with respective likelihood operators L^, 
and L'~^ . Operators for daughter nodes B, C are com- 
bined using the analog of pruning, the quantum prun- 
ing map p : Lin{H) (g) Lin{H) — > Lin{H) that pro- 
vides the parent operator = p{L^ ® L'^)i where 
pi = Tvb o AdUl^ o Edd o Ad{UB ® Uc). The map p 
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uses stochastic matrices M^{trc) = Ux oU*, depending 
on branch lengths for x — A,B, as given by the 
model employed, and the collective "diagonalizing map" 
£dd{-) = J2k Pk®Pk{-)Pk®Pk- Fig. m presents a quantum 
circuit realizing map By using its embedding ^r,r+i = 
^^®r-i^^^j^(gis-r jTQj. various values of r according to the 
topology of the binary tree, the pruning map /i is applied 
recursively to all cherries, and then higher up the tree. 
In this way we arrive at the tree likelihood operator L^^^ , 
which then is contracted with model's stationary density 
matrix — '^^T^iPi, to yield as a measurement result 
the site I likelihood L^) = Tr{LfJp'') = (Z'iJ , p'^'j . For 

the entire alignment, the tree (log) likelihood is (c.f. the 
identity Tr{AB)Tr{CD) = Tr{A ® C){B ® D) ) 

A 

L = maxlogTT ll^tl pA = maxlogTr{<i^tiW)pA), 

w I w 

1=1 

where pA = (p^)®^ is the product of A stationary density 
matrices. 

In fact this Heisenberg-like picture of updating the ob- 
servables (likelihoods), and finally contraction with the 
stationary density matrix to derive site and eventually 
alignment likelihoods, can be converted to a Schrodinger- 
like picture, using the observable-state duality, exempli- 
fied here by the trace cyclic property. Firstly note that 
the pruning map can be expressed as p.{L^ (8) L*-^) = 
iy^^£B{L'^), V = TrL^, where the positive stochas- 
tic map £b decomposes as £b = £Bpd ° AdUc with 
£Bpd{-) = ^ke^lk ^^Pk{')^ ^ probabilistic diagonaliz- 
ing map, with probabilities qj^ = v^^ {k\ UbL^U^ \k) . 
As the roles of and can be exchanged above 
with appropriate modification, {£b becomes £c etc), we 
note that /x is proportional to a stochastic map either 
way, and by duality it can be made to act on density 
matrices instead of likelihood operators. This is also 
true for embedded pruning maps pr,r+i, i-e. they will 
also be proportional to maps £B;r,r+i for the appropri- 
ate current likelihood etc. Then the tree likelihood 
operator L(^\ obtained by composing pruning maps, 
will eventually be described by pruning a final cherry, 
say with nodes Bf and C/, ie. Lj^^ = i'f~^£Bj{L'"^), 
Vf = TrL^f . Then the likelihood at site I is obtained 
as = ,yf-^Tr{£B,{LCf)pn - ^r^TriL^s E*s^^Xp^)), 
where the dual map £g of f acting on the density ma- 
trix is introduced. This situation is extended similarly to 
the likelihood of the entire alignment by assigning addi- 
tional site indices I to each likelihood operator, e.g. hf^ 
and l'^^ ^ as well as trace coefficients i/^ etc, to obtain 
L{T-w*) = max^lognti(^/)-'?^KiA('^ti^s,,)PA)- 

Here La = ®i=\ij^'i^') is the product of A different like- 
lihood operators, corresponding to final cherries of the 



respective trees, employed to construct tree likelihoods. 
Note that ®^^-fi*Bj , ^ collective factorized map that 
can be expressed in terms of a unitary dilation, and 
this would in principle be implemented by a Hamiltonian 
quantum model. 

In conclusion, this study lays the groundwork for sim- 
ulating, by quantum mechanical means, the probability 
tensors of multi-taxa systems, and for estimating the 
maximal likelihood of a phylogenetic alignment. With 
the tools developed here, prominent among problems for 
future investigations would be for example a quantum 
computational simulation of Steel's conjecturc^^ and its 
resolutions^. 
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